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Abstract 

In this paper, we present a fully discretized Calderon Calculus for the two dimensional 
Helmholtz equation. This full discretization can be understood as highly non-conforming 
Petrov-Galerkin methods, based on two staggered grids of mesh size h, Dirac delta distri- 
butions substituting acoustic charge densities and piecewise constant functions for approxi- 
mating acoustic dipole densities. The resulting numerical schemes from this calculus are all 
of order h 2 provided that the continuous equations are well posed. We finish by presenting 
some numerical experiments illustrating the performance of this discrete calculus. 

Key words: Calderon calculus, Boundary Element Methods, Dirac deltas distributions, Nystrom 
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1 Introduction 

In this paper we present a very simple and compatible Nystrom discretization of all boundary 
integral operators for the Helmholtz equation in a smooth parametrizable curve in the plane. 
The discretization uses a naif quadrature method for logarithmic integral equations, based on 
two staggered grids, and due to Jukka Saranen and Liisa Schroderus [13] (see also [15] and 
[2]). This is combined with an equally simple staggered grid discretization of the hypersingular 
operator, recently discovered in [8]. If the displaced grids used for the discretization of these two 
operators are mutually reversed, then it is possible to combine these two discretizations with 
a simple minded Nystrom method for the double layer operator and its adjoint. The complete 
set of operators is complemented with a fully discrete version of the single and double layer 
potentials. We will explain the construction of the discrete set and reinterpret it as a non- 
conforming Petrov- Galerkin discretization of the operators (using Dirac deltas and piecewise 
constant functions) to which we apply midpoint integration in every element integral. 

Once the semivariational form has been reached we will show inf-sup conditions for all 
discrete operators involved and consistency error estimates based on asymptotic expansions of 
the error in the style of [2j [5j [6] . We will finally state and sketch the proof of some convergence 
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error estimates. While some of the results, for individual equations (mainly based on indirect 
boundary integral formulations) had already appeared in previous papers, this is the first time 
that the entire Calderon Calculus is presented in its entirety. Let it be emphasized, that this is 
probably the simplest form of discretizing simultaneously all the potentials and integral operators 
for the Helmholtz equation in the plane and that the methods we obtain are of order two. Barring 
the conceptual difficulty of understanding the boundary integral operators, the methods have 
the simplicity of basic Finite Difference Methods and require no effort in their implementation: 
all discrete elements are described in full, natural data structures can be easily figured out from 
the way the geometry is sampled, and no additional discretization step (quadrature, assembly 
by element, mapping to a reference element) is required. The methods will be presented for 
the case of a single curve, but we will hint at its immediate extension to the case of multiple 
scatterers. 

In a final section devoted to numerical experiments, we will show how to use the methods 
for transmission problems and how to construct combined field integral representations. 



2 Calderon calculus for exterior Helmholtz boundary problems 

2.1 Potentials and operators 

Let r be a smooth simple closed curve given by a regular 1-periodic positively oriented parametriza- 
tion x = (x 1 ,x 2 ) : M -> T C M 2 . Let n(t) := (x' 2 (t), -x[(t)) be a non-normalized outward point- 
ing normal vector at x(t) G T. The domain exterior to T will be denoted Q + . As a reminder of 
the fact that we are taking limits from this exterior domain, the superscript + will be used in 
trace and normal derivative operators. 

Given 1-periodic complex- valued functions rj and tp, the (parametrized) single and double 
layer potentials are defined with the formulas 

(S,)(z) := jj 1 H<V(k\z- X (t)\)r,(t)dt, (1) 

fl{ 1) W«-x(t)|) (, -^yV (t)d t (2) 

for arbitrary z£l 2 \T. (Here is the Hankel function of the first kind and order n.) The 
single and double layer potentials define radiating solutions of the Helmholtz equation, namely, 
if XJ = Srj + D?/>, then 

AU + k 2 U = in fi+, W(z) • (Az) - ikU(z) = o(-4-), as |z| -> oo. (3) 

' ' vl z l 

Moreover, if U is a C 1 (0+) solution of © and we define 

^ = 7 +[/ := [[| r ox, \ = d+U:= ((W)|r°x)-n, (4) 

then [91 H5 

U(z) = (D <p)(z) - (S A)(z), z G + . (5) 

We note that the representation formula ([5]), depending on parametrized Cauchy data can 
be extended to any locally H 1 solution of ([3]). In this work we will restrict our attention to 
smooth solutions though. 



(DV-)W - f / 
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Associated to the layer potentials we have three integral operators. 



(Vr,)(s) :-- 
(K^)(s) :-- 



i 

4 



4 

ik 

T 



~ / ff{%|x( S )-x(i)|) 



^ (1) (fc|x( S )-x(t)|) 



(x(s) - x(t)) • n(t) 
|x(a)-x(t)| 

(x(t) - x(s)) • n(s) 
|x(a)-x(*)| 



n(t) dt, 



as well as the integrodifferential operator 

Wip := -(Yip' 

where 



(Vn^)(*) := 4 



x s 



x(t)|)(n(t)-n( S ))V(t)dt. 



(6a) 
(6b) 
(6c) 

(6d) 



The operators in ([6]) are respectively called single layer, double layer, adjoint double layer, and 
hypersingular operator. The operator W admits a different expression in terms of finite parts 
integrals (see |14| Lemma 2.5.6]), which is where its name comes from. 

Layer operators and potentials are related via the so-called jump relations [9| lllj. fTS] . namely, 
the exterior parametrized boundary values of the layer operators are given by the formulas 



7 + S 7] 
d+Srj 



-|r/ + J?7, d+Dip 



The matrix of operators 



I + K 

-W 



-V 



(7) 



(8) 



is the exterior Calderon projector. It follows from ([5]) and (J7J), that if (</?, A) are the parametrized 
Cauchy data (j3|) for a solution of ([3|), then C + (cp, A) T = (</?, A) T or, equivalently 







"il-K V 








"0" 


A_ 




W ±I + J_ 




A_ 








(9) 



Note that K and J are transposed of each other, while V and W are symmetric. 



2.2 Boundary integral equations for exterior problems 

We next summarize a collection of boundary integral equations leading to the solution of Q 
with a given boundary condition: 

7+C/ = /3 or dfcU = p 1 . (10) 

The data functions in the right-hand side of (|10p are 1-periodic functions and the boundary 
operators are those of ([!]). Recall that the Dirichet or Neumann exterior problem for the 
Helmholtz equation with Sommerfeld radiation condition at infinity are uniquely solvable. 

A direct method for solving the exterior Dirichlet problem starts in the representation formula 
©, equates (p = /3q, and then uses one of the two identities in ([9|) to set up an integral equation 
in order to find A. Similarly, for the Neumann problem, we impose A = /3\, and then use one of 
the equations in ([9]) in search of (p. The resulting integral equations are collected in Table [TJ 
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Dirichlet 


"\ T \ i TV 

VA = — ±cp + Kip 




= Po 


(dDOl) 












2 A + JA = —W/ip, 


<P 


= Po 




Neumann 


— \<p + Ky? = VA, 


A 


= Pi 


(dNOl) 












-Wp = ±A + JA, 


A 


= Pi 


(dN02) 



Table 1: BIEs for direct formulations. The representation formula is ([5]). All these equations 
are solvable. Uniqueness is discussed in Proposition 12.11 

An indirect method based on the single layer potential representation looks for U = Srj and 
then uses the expressions in the first column of ([7]) to set up an integral equation depending 
on which boundary data is known. Similarly, we can look for U = T)ip and use the boundary 
integral operators that appear in the right column of ([7]) to build an equation. These equations 
are gathered in Table [2j 



Dirichlet 


V?? = po, 




u = 


S 77 


(iDOl) 














+ = 


Po, 


u = 


DV 


(iD02) 


Neumann 


-\r\ + J77 = 


Pi, 


u = 


S 77 


(iNOl) 














wv = -Pi, 




u = 




(iN02) 



Table 2: BIEs for indirect formulations. The potential representation is given next to the 
boundary integral equation. Unique solvability of these equations is discussed in Proposition 
I2TT1 

Proposition 2.1 (See |12t Section 3.2]). Let Q be the domain interior to V . 

(a) Equations (dNOl), (iNOl), (dDOl), and (iDOl) are uniquely solvable if and only if —k 2 is 
not a Dirichlet eigenvalue of the Laplace operator in 0,. 

(b) Equations (dN02), (iN02), (dD02), and (iD02) are uniquely solvable if and only if —k 2 is 
not a Neumann eigenvalue of the Laplace operator in £1. 

The equations of Tables [T] and [5] involve the four operators of the matrix T> + in Q and 
their transposes. The operators in the first row of T> + are invertible when — k 2 is not an interior 
Dirichlet eigenvalue. The operators in the second row of T> + are invertible when —k 2 is not an 
interior Neumann eigenvalue. The precise Sobolev space setting where these equations are well 
posed will be explained in Section 14.11 In addition to these equations, the Calderon Calculus, 
given by the jump relations ([7]) and the identities (|9|) can be used to construct combined integral 
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equations and several other associated boundary integral equations, some of which are invertible 
for all values of k. 



3 The fully discrete calculus 
3.1 Matrix representation 

Let N be a positive integer, h := l/N, and let us consider the uniform grid in parametric space 

Si := (i — \)h, ti := ih, i € Z, 

thus defined so that U is the midpoint of the interval (sj, Sj+i). The following quantities will be 
all the geometric elements of T that will be used in the discrete Calculus: 



mi := x(ij), b { := x(sj), 

rij := hn(ti), ii := |n;| = h\x.'(U)\, = /i 2 x"(i;). 



(11) 



These quantities make up the main discretization grid. Note that they are defined for i 6 Z, 
modulo N. For practical reasons, we will need a discrete function n(i), that gives the next index 
in a rotating (modulo N) form, so that n(i) = i + 1 for i < N — 1 and n(N) = 1. We now take 
e £ (— l/2,l/2)\{0} and repeat the same construction with the displaced grid in parametric 
space: 

tf = (i + e)h, si = (i + e - \)h. 

The quantities mf, bf, nf, if, and sf are defined accordingly. They constitute the companion 
grid. 

Given column vectors r\ = . . . ,^at) t G C^, = (V>i, • • • , 4>n) T £ C N , we consider the 
discrete single and double layer potentials: 



N 



j'=i 

N 



D h (z)^ 



J^^z-m, 



j'=i 



in 



We also consider four N x N matrices V^, K^, J/j and W^, given by 

V,; 



^(fcK-mJI), 



K 



4tt£2 ' 



4 

sj • n; 



m, — m , 



4vr(£f)2' 

ik (i) e (m £ -mf) n 

'(fc|m? - m jD i m£ _ m£ | ' 



* = J, 

* = j, 



W 



where 



v n(i),n(j) + Vjj - V„ ( j)j - V ijn(j -) — A; (nf ■ n^V,^, 
Y l3 = l -H i i\k\hi-h ] \). 



(12a) 
(12b) 

(13a) 
(13b) 

(13c) 
(13d) 
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Note that the diagonal values in and Jh are defined using the limit values in the kernels of 
the integral operators K and J as \s — t\ — > 0. 

Remark 3.1. As can be seen from (|12p and (|13p . the structure of the matrices and operators 
does not remember where the discrete geometric data come from. The formulas (|12j) and Q13|) use 
discrete data {nij, nj, bj, £j, Sj} and {m|, n?, b|,0, s?}, sampled from the curve. It is immaterial 
whether these data have been sampled from a simple curve or several simple non-intersecting 
curves. The next-index function n{i) used in Wh has to be adapted to contain cycles of nodes 
showing the different connected components of the collection of curves. 

Discretization of the integral equations in Tables [1] and [2] is almost straightforward based 
on these matrices and potentials. The Dirichlet and Neumann data in (|10|) are discretized by 
vectors of samples: 

(3 :=((3 (t 1 ),...,p (t N )) T (ft (if),..., ft (**r)) T - ( 14 ) 

The different scaling of these vectors will be clear from the interpretation of these methods 
that we will give in Section 13,21 At this stage, it can be justified with some arguments of 
dimensional analysis, given the fact that j3\ corresponds to data of a derivative of the function. 
By the definition of the parametrized boundary operators (|1J), of the Cauchy data (fTUj) and of 
the discrete quantities (jlip . we can similarly write 

/3o = (U( mi ),...,U(m N )) T Pi = (VU(ui e 1 )-nl,...,VU(m e N )-n £ N ) T . 

The discrete direct methods use a representation formula 

U h (z) = D h (z)\-S h (z)<p (15) 

and one of the linear systems of Table [3l The discrete indirect methods appear collected in Table 
HI including the corresponding potential representation. 



Dirichlet 


Y h X = -\<p + K h tp 




= /3o 


(dDOlh) 












\X + J^A = -W h ip, 


<P 


= A) 


(dD02h) 


Neumann 


-±(p + K h ip = \ h X, 


A 


= A 


(dNOlh) 












-W h ip = ^\ + JhA, 


A 


= Pi 


(dN02h) 



Table 3: Discrete direct methods, with representation formula (fT5"jh 

3.2 Reinterpretation as non-conforming Petrov-Galerkin methods 

Our method can be understood as a collection of non-conforming Petrov-Galerkin methods with 
a very simple quadrature rule for approximating any integral appearing in the scheme. The 
basic idea is the following: the input of D (and therefore W and K) will be approximated with 
a piecewise constant function on the main grid; the input of S (and therefore V and J) will be 
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Dirichlet 


VhV = /3q, 






= ^>hV 


/■r\ni 1 \ 

(lDOlh) 














1 _ /, 1 TV" _ /. 

2^ + <^hW = 


Po, 


T J 


= ^>hW 


(lUUznj 


Neumann 


+ = 


/3i, 




= ShV 


(iNOlh) 














W fe -0 = -A, 




u h 




(iN02h) 



Table 4: Discrete indirect methods. 



approximated with a linear combination of Dirac deltas on the companion grid; tests related to 
Dirichlet problems will be carried out by Dirac deltas on the main grid; test related to Neumann 
problems will be done with piecewise constants on the companion grid; finally, all integrals will 
be broken into subintervals of the grid and approximated with a midpoint rule. 

In order to write the methods of Section 13.11 in the form where we will develop their con- 
vergence analysis, we need to define some new discrete elements. First of all, we consider the 
(periodic) Dirac delta distribution S z at a point z. Its action on any periodic function that is 
continuous around z will be denoted {S z ,p} = {p, 5 Z } := p(z). Given an interval I C R, we will 
denote by \i the periodized characteristic function of /, i.e., the characteristic function of the 
set I + Z. We then consider four discrete spaces 

S h := span{x( Si _ 1)Si ) : i = l,...,N}, S h)E := span{x( s |_ 1)S |) : i = l,...,N}, 
1 := span{(5 ti : i = 1, . . . , N}, S^ £ := span{<5^ : i = 1, . . . , N}. 

For elements of these spaces we will identify the vector of their coefficients -with respect to 
the basis that has been used to define the space-, using the same letter in boldface font. For 
example, 

N 

s^ 1 3 ph = y y fijStj < — > n = (pi, ■ ■ -,pn) t e c^. 

The two discrete operators 

N N 

Qh'p ■= h £ P(h)5 t] Ql] e p := h £ pflfy 
i=i i=i 

complete the collection of elements needed for a more variational description of the discrete 
Calderon Calculus. They will be used to denote midpoint quadrature approximations. For 
example, 

N i 

{Q^p,(f>} = hJ2p(tj)<f>(t3)~ / P(t)<f>(t)dt. 

3=1 J ° 

The discrete potentials (fl~2j) can be easily described in this language: 

Sj£ 3i)i4 S h {-)r] = Sr} h , S h 3ip h ^ D h (-)ijj = DQ~V/i- 

Observe how in the double layer potential we are just applying the midpoint rule to approximate 
Diph, while no additional integration is needed in the already fully discrete expression for Srjh- 
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The matrices (|13p have their variational counterparts as bilinear forms: 



Sh 1 x S h*e B (vhiVh) 1 — > y(fj-h,Vh) ■= {vhiVvh} = m t V/j»7, 

S^ 1 x S h 3 (fJ-h,4>h) 1 — >■ k(fj, h ,iph) ■= {fJ>h,KQh 1 iph} = M T Kfci/'> 

Sh,eX S^3 ((f>h,T]h) ' ► }(4>h,T]h) '■= {Qh^h, «%} = & T Jhr), 

S h>£ x S h 3 ((j) h ,ip h ) — ► ™(<l>h,iph) ■= T WkV>- 
The bilinear form w can be understood as follows 

w(^,^) = V<} - k 2 {Qll<t> h , VnQ^V/,}, 

just by noticing that x'( s . x s ) = — an d that a change of sign has to be applied to the 
leading integrodifferential part of W (see (|6dj) ) when changing the differentiation to the test 
function. The rationale behind this choice of spaces can be observed in the matrix of operators 
T> + in ([9|). As trial spaces we are considering 5^ x Sj^ , while the rows of T> + are respectively 
tested with S^ 1 and Sh, £ - This means that the operators of the second kind (±^I + K and 
=t^I + J) are discretized on a single grid (each of them on a different grid though), while the 
operators of the first kind (V and W) use two grids. This is actually a requirement due to the 
fact that the kernels of V and V n cannot be evaluated in the diagonal s = t, where they have a 
logarithmic singularity. Once this choice of trial and test spaces has been taken as a first step in 
the discretization of the four operators in ([9|), midpoint integration is applied to all remaining 
integrals. The operators Qj^ 1 and are used as a way of enforcing full discretization of every 
operator acting on a piecewise constant function. 

To describe variationally the equations in Tables [3] and H] we first cast the data function (j3o 
for the Dirichlet problem and /3\ for the Neumann problem) in the discrete spaces 



N N 
3=1 3=1 

so that their coefficients coincide with the sample vectors Q14|). The equations (dNOlh) corre- 
spond then to writing Aft, = f3±, solving 

ifh^S h s.t. -l{lih><Ph} + Hl J 'h,<Ph)=v(iihAh) Vft e S,; 1 , (16) 

and finally using Uh = DQ^tph — SA^ as discrete representation formula. The indirect method 
(iN02h) corresponds to solving 

ijlh^Sh S.t w(0 ft ,^fc) = -{Pi,<f>h} = -{Qh^hiPl} V4>h£S h , £ , 

for a potential representation Uh = DQ^t/V The indirect method (iDOlh) is equivalent to 
solving 

T] h eS^ e s.t. v(/i ft ,77 h ) = {/J,h,Po} = {VhiPo} V/ift 6 S^ 1 . 

The remaining five discrete equations in Tables and U] can be easily rewritten using these same 
elements. 
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4 Numerical analysis 



4.1 Stability 

Analysis of the methods in Section 13.11 is carried out in the form given in Section 13. 2\ in the 
frame of periodic Sobolev spaces. For s£lwe define the space H s as the completion of the 
space of trigonometric polynomials span {exp(2-7rim • ) : m G Z} with respect to the norm 

||p||2 = |^(0)| 2 + V |m| 2s |p(m)| 2 , p(m):= [ p(t) exp(-2vrzmt) dt. 
, n Jo 

m^O 

An extensive treatment of these spaces can be found in [2] . The operators ([6]) can be extended 
to act on all Sobolev spaces H s . In particular, the following result holds (see Table 2.1.1] 
and [121 Section 3.2]). 

Proposition 4.1. The operators 

± \ + K, ±± + J : H s -> H s , V:H S ^H S+1 , W : H s -»• fl* -1 (17) 

are bounded for all s. If, in addition, —k 2 is neither a Dirichlet nor a Neumann eigenvalue of 
the Laplacian in £1 ( cf. Proposition \2.1\) . then all of them are invertible. 

Proposition 4.2. Assume that —k 2 is neither a Dirichlet nor a Neumann eigenvalue of the 
Laplacian in Q, and let e G (—1/2, l/2)\{0}. Then there exist positive numbers cy, ck, cj, cw > 
so that for all h small enough 

mi sup j — j: — - — > cy, (18) 

\±^{p: h ,tp h } + k(n h ,tp h )\ 
mf sup ^— n — n — — — n > c K , (19) 

o^Vhe5 ho ^ g5 -i llM/ill-illwllo 

| ± +j(^fc,T/ft)| 

mf sup — n > cj, (20) 

o^r,h£S^ s o/^es feiS ||^ ||o ||% II -l 

. , Iw^h,^)! , . 

mi sup n t n n , n > cw- (21) 

o^ h eS h o^4, h es h}£ H^/iHoll^/illo 

The constants can depend on e. 

Proof. Condition (|18p was proved in [21 Proposition 8], although it is based on a stability result 
(phrased in different terms) given in [13] . Condition (|2ip has been proven in [SJ Theorem 1]. 
With minor modifications, the proof of Theorem 2] can be used to prove (|20p . It is then 
easy to note that this result would also hold for the spaces S*^ 1 and Sh (it all amounts to 
displacing the grid for both test and trial functions). Then, by an easy transposition argument, 
(USD holds. □ 



The value e = is not a practicable option for the choice of the grids: in this case both 
grids coincide and we are obliged to evaluate the singular kernels in their diagonal. The choices 
e = ±1/2 lead to a discretization of V (they give the same one) that is not stable, i.e., the 
inf-sup condition does not hold. The proof of the inf-sup condition for the discretization of W 
in [8] requires also that e ^ ±1/2, because it is based on the result for V, although numerical 
evidence points to this being just a technical restriction, which is not in the case of V. Note 
finally that dependence of the methods on e is 1-periodic. 
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4.2 Consistency analysis via asymptotic expansions 

We next study the consistency of the approximation of the bilinear forms associated to the 
four operators ([6]) by their discrete counterparts, as well as the approximation of the identity 
operators that appear in the equations of Tables [1] and [21 The consistency error analysis is 
carried out by comparison with a quasioptimal projection of the corresponding unknown (the 
input of the integral operator) in the discrete space. These projections are defined by matching 
the central Fourier coefficients: 

S hi 9 D hiv, D^Jl(m)=v(m), -N/2<m<N/2, 
S h 3 D h il), 5J(m)=^(m), -N/2 < m < N/2. 

The operator Dh was studied in [I], while proceeds from [2]. It is proved in those references 
that 

\\D h fj; - f/;\\ a < C s>r h r - s U\\ r s < r < 1, s < 1/2, (22a) 
\\D^ £ r) - rj\\ s < Ch r - S \\r]\\ r s < r < 0, s < -1/2. (22b) 

Proposition 4.3. For all n E H 3 and tp E H A it holds 

\{<f>h,D^n}-{</> h ,Q^rj}\ < Ch 3 \\ V \\ 3 \\<t> h \\o, V0 fc E S h , E , 
\{fi h ,D h i>}-{v h ,i>} + h 2 ±{ii h ,iJ;"}\ < CT^IMUKH-i, V/i ft eS h \ 

The constants in the bounds are independent of e. 



Proof. The second expansion follows from [3 Theorem 7]. To prove the first one, note that by 
[21 Lemma 5] 

L) h^ r l ~ Qh^e 1 ! = Qh}e E hi], where E h n := rj(m) ex.p(2mm-) - n. (23) 

-f< m <f 

A direct computation (see also [U Lemma 9]) shows then that 

\{4>h,QhlE hV }\ < 11^11011^110+ {<t>h,QhlE h ri}- [ 4> h (t)(E h v)(t)dt 



o 

< ll^Hodl^Ho + nhWE^h) < C^H^HolMla 

where the last inequality follows from the fact that H-EX^Hs < C/i* _s ||?7||j for all t > s [L 
Theorem 8.2.1]. □ 

For simplicity, in what follows we will write P E £(n) when P is a periodic pseudodifferential 
operator of order n, i.e., P : H s — > H s ~ n is bounded for all s. 

Proposition 4.4. There exists Pk E £(1) so that for all n E H 3 and ip E H 4 , 

Wh,D^ri)-{Q^</> h ,3 V }\ < C^Hr/Hall^llo, V0 fc E S h , £ , 
Un h ,D h ^)-{^Kii}-h 2 {^ h ,V^}\ < C^H^mil/ifcll-i, V^E^" 1 . 

The coefficient P^ and the constants in the bounds do not depend on e. 

Proof. We refer to [7], where similar expansions are derived. □ 
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The study of the approximation properties of V and W is strongly influenced by the param- 
eter e. We write 

Ci(e) := ^log(sin 2 (7re)) C 2 (e) := \ J* d(t)dt, 
and note that Ci(±l/6) = 0. 

Proposition 4.5. There exists a smooth function a v and operators L v G £ (1), G £(2), 
L>, L 2 , 6 G 5(3) such that for all r] e H 3 and ifr G 

K/M)- /^) -{WmVt?} - /iCi(e){/i^,a v r?} 

-/i 2 C2(e){w t ,L v 7?}| < C^II^Hsll^ll-i, V/i fc G S^ 1 , 
|w(^,D ft ^) - {Qj^Mty} - ^i(e){Qj^ & ,L^} 

-h 2 {Qh&h, (C2(e)\l a + L 2b )V}| < C^HVlUll^llo, V<^ G 5 fc , e . 

Proof. The first expansion is given in [21 Theorem 7], while the second one is proved in [8j 
Proposition A. 4]. □ 

The key fact at this point is that by letting e = ±1/6 all the expansions start at h 2 . This will 
be crucial since, as we will see in the next subsection, we can identify the order of the method 
with the first power of h appearing in the consistency expansion. The relevance of identifying 
the h 2 term of the asymptotic expansion of the consistency error in Propositions 14.31 14.41 and 
14.51 is related to the possibility of moving from the norms given by the inf-sup conditions in 
Proposition 14.21 to stronger norms when producing estimates of the convergence error. (See 
Theorem 14.71 below.) 

Remark 4.6. If in Propositions \4-4\ an d \4-5\ we on ^V assume that r\ G H 2 and -0 G H 3 , 

and we eliminate the h 2 term from the left-hand side of the bounds, then the result holds with a 
bound of the form Ch 2 \\rj\\2 or Ch H^lb- 

4.3 Convergence estimates 

We collect in this subsection the convergence results for the all numerical schemes presented in 
this paper. 

Theorem 4.7. Assume that k satisfies the hypothesis of Proposition \4~^ and e = ±1/6. Let 
((ph,Xh) £ SZ* x Sh be the pair associated to the solution (A, cp) of any of (dDOlh), (dD02h), 
(dNOlh) or (dN02h). Then 

\\<p h - D h ip\\ + \\X h - A||_i < C7i 2 (|M| 3 + ||A|| 2 ). 

Moreover, 

max\(pj - A)(i/)| + max\h~ 1 X j - /h(tj)\ < Ch 2 {\\<p\\ A + ||A|| 4 ). 
j j 

Proof. We will only show the case (dNOlh), all others being very similar. Using the variational 
representation of (dNOlh) in (|16p . we can write 

-j{fj, h ,ip h } +k((j, h ,ip h ) = {nh, ~l<p + Kip} + {fj, h ,V(Q^ £ \ - A)} Mp, h G S^ 1 . 
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Using now Propositions 14.31 (second bound). 14.41 f second bound) and 14.51 (first bound) -see also 
Remark 14.61 - and (|23p . it follows that 

I - ^{fJ-h,<Ph ~ D hV} + H^h,fh - D h tp)\ < C7/i 2 1 1 /x^ 1 1 _ a_ ( 1 1 1 1 3 + ||A|| 2 ) 

+M»h,QhiE h \)\ V/ihG^ 1 . 

Using [2, Lemma 13] and the fact that 

\\Qh*v\\-i<CQ\v\\o + h\\r,\\x), 

(see [21 Lemma 6]) we can prove that 

\v^ h ,Q^ E E h X)\ < C||// h ||_i(p fc A||o + HEhXWt) < C^llMftH-illAHa fyh G S^ 1 . 

Therefore, by Proposition 14.21 the bound for ||^ — Dh^pWo follows. The bound for 

||A - AfcH-! = ||A - Qfc\\-i < || A - D^A||_! + \\D& - Q£\\\-i 

follows from ()22bj) and (j23|) . The uniform estimates require including the h 2 term of the consis- 
tency error expansion: see [21 Corollary 11] and [51 Theorem 6.4] for very similar arguments. □ 



Theorem 4.8. Assume that k satisfies the hypothesis of Proposition and e = ±1/6. Let 
iph G Sh be associated to the solution ip of (iD02h) or (iN02h) and let rjh G be associated to 
the solution r) of (iDOlh) or (iNOlh). Then 

\\D h i/> - Vhllo < Ch 2 \\t/j\\ 3 \\ V - TlhW-! < Ch 2 \\ V \\ 2 . 

Proof. The proof is very similar to the one of Theorem I4.71 The absence of integral operators 
in the right hand side makes the arguments slightly simpler. □ 

In all cases it is possible to prove that the estimates can be transferred to the computation of 
potential, with the direct representation (|15|) in the case of direct method, or the associated layer 
potential in the case of indirect methods. In all cases, we can prove \U(z) — Uh(z)\ < C(z)h 2 . 

Remark 4.9. If we take e ^ ±1/6, the methods involving or are of order one. 



5 Experiments 

In the following experiments we consider a single elliptical obstacle with boundary 

1 ( X _ 0.1)2 + {y _ 0.2)2 = 1. 

We will check solutions in two observation points inside the ellipse xi = (0.2,0.4) and X2 = 
(—0.2, —0.4). The examples will use more complicated integral equations than those explained 
in the previous sections, in order to put the discrete Calderon Calculus to a more demanding 
test. 



12 



5.1 A transmission problem 

Consider the coupling of the exterior Helmholtz equation ([3]) with an interior equation with 
different wave number 

AV + {k/c) 2 V = in O 
(here c > 0) and transmission conditions 

1 +U + f3 = 7 -V, d+U + P^ad-V 

(with a > 0). Data are taken so that the exact solution is 

U(z) = H^\k\z-xo\), V(z) = exp(i(k/c)z • d) x € ft, |d| = 1. 

We use the symmetric formulation of Martin Costabel and Ernst Stephan [1] (see also |1U]). The 
main unknowns are (p~ = "f^V and A~ = ad~V. The system they satisfy is 



W k + aW k/c J fc + J fc/c 








W fe il + J fc " 




" A) " 


-K k - K k/c V fc + ^V fc/c _ 




_ A" _ 




. §1 - Kfc Vfc 




. ft . 



(24) 



where we have tagged the integral operators with the corresponding wave number. The potential 
representation for the interior and exterior fields is 

U = -$ k (\- - fc) +T> k (<p- - fa), V = a- 1 S k/c \- -D k/C <p-. (25) 

Discretization is carried out by simply substituting the elements of (|24|) and (|25|) by their discrete 
counterparts: the data functions are sampled with (|14p . the integral operators are build with 
(|13|) and the potentials with (|12|) . We solve and tabulate the following errors: 

E£ := max iTr^ - ad+V(fi)\ E£ := max\<pj - 7~V(t,-)l 

EX :=max|^(x,)-V(x,)| 

These experiments are reported in Tables and The parameters are k = 3, c = 2/3 and 
a = 3/2. 



N 


error 


e.c.r 


N 


error 




e.c.r 


10 


4.6842£(+000) 




10 


5.8671£(- 


-001) 




20 


1.2470£7(+000) 


1.9093 


20 


1.9970E7(- 


-001) 


1.5542 


40 


3.7207£(-001) 


1.7448 


40 


4.9104£(- 


-002) 


2.0246 


80 


9.4663£(-002) 


1.9747 


80 


1.2376£(- 


-002) 


1.9883 


160 


2.3768£(-002) 


1.9938 


160 


3.1081£(- 


-003) 


1.9934 


320 


5.9518£(-003) 


1.9976 


320 


7.7699£(- 


-004) 


2.0001 


640 


1.4886£(-003) 


1.9994 


640 


1.9423E(- 


-004) 


2.0001 



Table 5: Errors E^ (left column) and E^ (right column) for the Transmission Problem in Ex- 
periment 1. 
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N 


error 




e.c.r 


10 


1.8729£(- 


-001) 




20 


2.0779£(- 


-002) 


3.1721 


40 


4.0885£(- 


-003) 


2.3455 


80 


9.6559£(- 


-004) 


2.0821 


160 


2.4527£(- 


-004) 


1.9770 


320 


6.1837£(- 


-005) 


1.9878 


640 


1.5527£(- 


-005) 


1.9937 



Table 6: Error EY (potential solution V at two interior observation points) for the Transmission 
Problem in Experiment 1. 

5.2 Burton-Miller integral equation 

Consider now the exterior Helmholtz equation ([3]) with boundary condition + 7 kmc = 0, 
where AC/i nc + k 2 U- mc = in a neighborhood of the interior domain Q. The well known Burton- 
Miller integral equation Section 3.9] is 

\i + ^ + cVi = dJJ- mc + c 1 U inc . (26) 

The exterior normal derivative can be computed after solving this equation and there are two 
potential representations of the solution 

A = i - d n U inc U = -S£ = -SA - D 7 C/ inc . (27) 

The value c = — ik is the usual choice in (|26p . For this value, the equation ()26p is uniquely 
solvable independently of the frequency. Since S£ = U- mc in the interior domain, we compare 
errors 

E^ := max jS/^x^ - U inc (x e )\ 

We also compare the density £ with the solution of Problem (dDOlh) (Table ED computing the 
compared error 

Ej :=max| h^Xj - (hT^ - d n U inc (tj)) \- 

(dDOlh) Burton-Miller I.E. 

In our numerical experiments we have taken ?7 mc (x) = exp(ikd ■ x), i.e. an acoustic plane wave, 
with direction given by the unit vector d = (1, l)/y/2 and wave number k = 2. The results are 
gathered in Table [7J 

5.3 Conclusions 

We have presented a collection of compatible discretizations of the two potentials and four 
boundary integral operators associated to the Helmholtz equation on smooth parametrizable 
curves in the plane. We have shown discrete stability of the discrete versions for all the operators 
in absence of resonances. We have also given convergence estimates for eight integral equations 
that solve the exterior Dirichlet and Neumann problems, with direct and indirect boundary 
integral equations. Finally, we have tested the methods in more complicated cases, such as 
systems of boundary integral equations arising from transmission problems and combined field 
integral equations. 
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N 


error 




e.c.r 


N 


error 


e.c.r 


10 


1.7205£(- 


-001) 




10 


7.6790£(+000) 




20 


3.6082£(- 


-002) 


2.2535 


20 


1.8790£(+000) 


2.0310 


40 


1.1990£(- 


-002) 


1.5894 


40 


4.1656£(-001) 


2.1734 


80 


3.7936£(- 


-003) 


1.6602 


80 


8.5219£(-002) 


2.2893 


160 


1.057LE(- 


-003) 


1.8435 


160 


1.4703E(-002) 


2.5351 


320 


2.7581£(- 


-004) 


1.9384 


320 


2.2452£(-003) 


2.7112 


640 


7.2185£(- 


-005) 


1.9339 


640 


7.1749£(-004) 


1.6458 



Table 7: Errors (left columns) and E? (right columns) for the Burton-Miller integral equation 
in Experiment 2. 

References 

[I] Arnold, D.N., A spline-trigonometric Galerkin method and an exponentially convergent 
boundary integral method, Math. Comp., 41 (1983), 383-397. 

[2] Celorrio, R., Dommguez, V., and Sayas, F.-J., Periodic Dirac delta distributions in the 
boundary element method, Adv. Comput. Math., 17 (2002), 211-236. 

[3] Colton, D. and Kress, R. Integral equation methods in scattering theory. John Wiley & Sons, 
New York, 1983. 

[4] Costabel, M. and Stephan, E.P., A direct boundary integral equation method for transmis- 
sion problems, J. Math. Anal. Appi, 106 (1985) 367-413. 

[5] Dommguez, V. and Sayas, F.-J., Local expansions of periodic spline interpolation with some 
applications, Math. Nachr., 227 (2001), 43-62. 

[6] Dommguez, V. and Sayas, F.-J., Full asymptotics of spline Petrov-Galerkin methods for 
some periodic pseudodifferential equations, Adv. Comput. Math., 14 (2001), 75-101. 

[7] Dommguez, V., Rapun M.-L., and Sayas, F.-J., Dirac delta methods for Helmholtz trans- 
mission problems, Adv. Comput. Math., 28 (2008), 119-139. 

[8] Dommguez, V., Lu, S., and Sayas, F.-J., A Nystrom method for the two dimensional 
Helmholtz hypersingular equation, Submitted. arXiv: 1210.4582 

[9] Hsiao, G.C. and Wendland, W.L., Boundary Integral Equations, Applied Mathematical Sci- 
ences. Springer, 2010. 

[10] Laliena, A.R., Rapun, M.-L., and Sayas, F.-J., Symmetric boundary integral formulations 
for Helmholtz transmission problems, Appl. Numer. Math. 59 (2009), no. 11, 2814-2823. 

[II] McLean, W., Strongly elliptic systems and boundary integral equations. Cambridge Univer- 
sity Press, Cambridge, 2000. 

[12] Nedelec, J.C., Acoustic and Electromagnetic Equations: Integral Representations for Har- 
monic Problems. Vol. 144 in Applied Mathematical Sciences. Springer, 2001. 

[13] Saranen, J. and Schroderus, L., Quadrature methods for strongly elliptic equations of neg- 
ative order on smooth closed curves, SIAM J. Numer. Anal., 30 (1993), 1769-1795. 



15 



[14] Saranen, J. and Vainikko, G., Periodic integral and pseudodifferential equations with nu- 
merical approximation, Springer, 2002. 

[15] Sloan, I.H. and Burn, B. J., An unconventional quadrature method for logarithmic-kernel 
integral equations on closed curves, J. Integral Equations Appl, 4 (1992) 117-151. 



16 



